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A  two-dimensional  viscoplastic  beam  element  was  developed 
for  the  computer  program  FEAP.  The  element  is  suitable  for 
small  displacements  or  moderate  rotations  with  negligible  shear. 
The  element  is  formulated  in  moment-»thrust  space  which  will 
improve  modeling  of  concrete  beams.  The  viscoplastic  formula¬ 
tion  provides  time  dependent  effects  such  as  the  influence  of  th< 

strain  rate  unnn  tho  violH  ctronnth  % _ 


1.  Introduction 

This  report  describes  a  two  dimensional  viscoplastic  plane  frame  beam  element  written  for 
the  computer  program  FEAP  [1],  and  is  a  summary  and  extension  of  work  presented  by  the 
first  author  in  [2].  The  element  is  suitable  for  the  analysis  of  structural  beam  elements  with 
rectangular  or  I  shaped  cross  sections  undergoing  small  displacement  or  moderate  rotation 
deformations  for  which  the  influence  of  shear  strain  is  negligible.  The  effects  of  both  static  and 
dynamic  loadings  may  be  considered. 

The  viscoplastic  constitutive  assumption  is  by  definition  a  time  dependent  plasticity  for* 
mutation,  in  which  the  rate  of  straining  affects  the  ultimate  yield  strength.  This  is  of  particular 
interest  in  the  analysis  of  structures  subjected  to  high  rate  blast  loadings.  It  is  also  possible  to 
recover  inviscid  plasticity  through  a  penalty  approach  by  appropriate  definition  of  the  material 
properties. 

1.1.  Report  Layout 

Section  Two  describes  the  development  of  the  weak  form  of  the  equations  of  equilibrium 
for  the  cases  of  small  displacement  and  moderate  rotation  deformation.  The  constitutive  theory 
and  its  adaptation  to  a  beam  formulation  is  presented  in  Section  Three.  Interpolations  for  stress 
resultants  and  centroidal  axis  displacements  are  introduced  and  used  to  assemble  the  governing 
algebraic  equations  for  the  element  in  Section  Four.  Numerical  examples  are  presented  in  Sec¬ 
tion  Five.  Appendices  I  and  II  contain  explicit  descriptions  of  the  component  element  matrices 
and  of  the  Newmark  time  integration  scheme  used  for  solution  of  the  nonlinear  dynamics  prob¬ 
lem.  A  description  of  the  element  input  data  is  found  in  Appendix  III. 
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2.  Week  Fora  Of  The  Eqeatiens  Of  EialllMm 

In  accordance  with  the  Bernoulli  -  Euler  hypothesis  that  plane  sections  remain  plane  and 
do  not  route  with  respect  to  the  centroidal  axis,  it  is  possible  to  define  the  state  of  stress  in  a 
beam  using  three  stress  resultants  --  the  axial  force  p,  shear  force  v  and  bending  moment  m. 

For  a  beam  subjected  to  concentrated  loads  applied  quasi-statically  at  the  ends  only,  the 
momentum  balance  principles  require  that  net  force  and  moment  on  a  differential  segment 
must  vanish.  The  positive  conventions  are  taken  as  in  Figure  1.  In  the  case  of  small  displace¬ 
ments,  the  equations  for  axial,  transverse  and  roUtional  equilibrium  take  the  form 

P,„  -  0  (la) 

v,x  -  0  (lb) 

mtM  “  v  (lc) 

where  (  ),,  indicates  differentiation  with  respect  to  variable  x. 

As  no  specific  allowance  is  made  for  the  effect  of  deformation  due  to  shear,  eqs.  (lb)  and 
(lc)  may  be  combined  to  yield  the  following  reduced  set  of  equilibrium  equations. 

P„  -  0  (2a) 

-  0  (2b) 

These  relations  must  be  identically  zero  within  the  domain  ft  *  of  each  element,  and  thus 

the  weighted  integral  form  of  eqs.  (2)  must  also  vanish  for  arbitrary  weight  function  W'.  The 

bold  type  indicates  a  vectorial  or  matrix  quantity,  and  superscript  /  indicates  the  transpose 
operation. 

/"'IXj'O'-o  «> 


Letting  u  represent  the  axial  displacement  and  w  the  transverse  displacement  in  a  given 
beam  element  of  length  I,  the  weight  function  W'  may  be  chosep  so  (^aL,V)rUt^l^ork  quanti- 

fie*  result  «  *  '  **  I 


ties  result.  ^  ^  ^ 

W'  =  Uu«lw)  (4) 

1  i\  .t/t  »■*!  j 

Substitution  into  eq.(3)  followed  by  integration  by  iparts  results  th^Gplqi&in.  form  of  the 
equations  of  equilibrium.  * 

. j 

^  {  8u,,  8w,„  )  j^jdx  -  '!  (5) 

The  boundary  terms  in  this  equation  form  the  jnatural  boundary  conditions  for  the  prob¬ 
lem,  and  hence  may  be  replaced  by  the  inner  product  of  the  virtual  nodal  displacements  81! 
with  the  applied  nodal  loads  q.  if 

It  will  be  necessary  to  linearize  this  equation  in  order  that  an  iterative  Newton-Raphson 
solution  scheme  may  be  used.  In  anticipation  of  the  rate  dependence  in  the  constitutive  theory, 
a  time  stepping  form  of  linearization  is  Employed.  Subscripts  n  represent  the  time  step 
number,  superscripts  i  the  iteration  counter  within  a  step,  and  A  an  incremental  quantity.  Fol¬ 
lowing  the  procedure  employed  in  [2],  the  linearized  form  of  the  small  displacement  equili¬ 
brium  equations  is  found  to  be 


•u'q.+t 


(6) 
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For  the  case  of  moderate  rotations,  in  which  the  transverse  deflection  w  is  of  the  order  of 
the  beam  depth,  the  equilibrium  equations  are  written  [3] 


* 

« 

© 

(7a) 

V,x  “  (P*,x),M 

(7b) 

«,x-  V 

(7c) 

As  before,  eqs.  (7b)  and  (7c>rgre  combined,  yielding 

P.jr-0 

(8a) 

m,„  -ipw,x),x  -  0 

(8b) 

The  Galerkin  form  for  the  moderate  rotation  equations  of  equilibrium  is  then  obtained  as 

in  the  previous  case. 

|(»K,x  +  *W.xO 

-  8up|f  +  —0 

(9) 

Linearizing, 

||(8u,x  +  8w,xw,x)  8w)XX|{  '  +  Am}„+1 

+  /|8w,xpAw,Jtjdr-  bU'q^i 

(10) 

The  second  integrand  in  eq.(10)  contains  a  nonlinear  term  involving  the  product  of  axial 
force  p  with  slope  w,x,  which  is  known  as  the  geometric  stiffness ,  Kc 


3.  Constitutive  Theory 

The  viscoplastic  constitutive  theory  [4]  postulates  that  the  total  strain  rate  may  be  addi- 
tively  decomposed  into  the  sum  of  elastic  and  inelastic  components, 

«-€*  +  «'  (11) 

with  the  elastic  strain  rate  being  related  to  the  stress  rate  by  the  elastic  compliance. 

-  C*<r  (12) 

A  closed  yield  surface  /  is  defined  in  stress  (or  stress  resultant)  space,  such  that  for  states 
of  stress  lying  within  the  surface,  the  material  response  is  purely  elastic.  The  inelastic  strain 
rate  is  then  given  by 

«'.y<*(/-)>!£  (13) 

where 

y  fluidity  parameter 

4(f)  function  chosen  to  represent  the  particular  material 

<4(f)>  operator  such  that  <4(f)>  -  0,  /  <  0 

<4(f)>-4(fh  />  0 

It  is  apparent  that  the  inelastic  strain  rate  vector  is  directed  parallel  to  the  gradient  of  yield 
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surface  /,  in  accordance  with  Drucker‘s  stability  postuiate(S].  In  contrast  to  the  classical  plasti¬ 
city  theory,  stress  states  outside  of  the  yield  surface  are  allowed,  but  only  with  non-zero  rates  of 
strain.  This  will  result  in  a  viscoplastic  flow  phenomenon.  Similarly,  if  the  displacement  field  is 
fixed  by  external  constraints,  the  stresses  will  relax  with  time  onto  the  yield  surface,  provided 
an  equilibrium  solution  exists  there. 

The  rate  of  viscoplastic  flow  or  stress  relaxation  depends  upon  the  stresses  indirectly 
through  the  function  4(/),  which  for  the  current  application  is  taken  to  be  a  power  law[4]. 

*(/)  =  /—  (14) 

The  magnitude  of  inelastic  strain  rate  is  also  affected  by  fluidity  factor  y.  For  a  given 
state  of  stress,  if  y  is  doubled,  the  inelastic  strain  rate  and  hence  inelastic  strain  increment  are 
also  doubled.  By  setting  y  sufficiently  large,  it  is  possible  to  model  invisdd  plasticity.  This 
forms  the  penalty  approach  alluded  to  in  the  introduction.  The  user  is  warned  that  y  may  not 
be  increased  without  bound,  as  indefinite  stiffness  matrices  will  result.  Selection  of  reasonable 
values  for  y  is  illustrated  in  the  section  on  numerical  results. 

As  a  first  step  in  adapting  the  constitutive  theory  to  the  beam  element,  it  is  necessary  to 
define  the  yield  surface  /.  Consistent  with  the  restrictions  stated  in  the  introduction,  it  is 
sufficient  to  consider  only  the  axial  normal  stress  in  determining  /.  The  cross  section  is  con¬ 
sidered  to  be  fully  plastified  when  all  fibers  have  yielded  in  tension  or  compression.  Intermedi¬ 
ate  states  are  not  included. 

The  stress  distribution  consists  of  two  rectangular  blocks,  one  in  tension,  the  other  in 
compression,  and  is  parameterized  by  a,  the  distance  from  the  centroidal  axis  of  the  section  to 
the  neutral  axis  of  the  stress  distribution.  (See  Figure  2)  When  a  —  0,  there  is  no  axial  force 
acting  on  the  section,  and  when  a  -  ±  A/2  there  is  no  bending  moment.  These  represent  the 
extrema  of  the  axial  force  -  bending  moment  interaction  curve,  or  yield  surface. 

It  is  necessary  to  consider  two  separate  cases,  one  in  which  the  neutral  axis  lies  in  the 
web,  and  the  other  in  which  it  lies  in  one  of  the  flanges.  This  is  facilitated  [7]  by  splitting  the 
stress  blocks  into  symmetric  and  anti-symmetric  parts  with  respect  to  the  centroidal  axis  of  the 
section.  The  symmetric  part  contributes  only  to  axial  force,  and  the  anti-symmetric  part  only  to 
bending  moment. 


Case  1 :  The  neutral  axis  in  web.  (a  <  A/2  -  t/) 

The  stress  resultants  p  and  m  are  first  expressed  in  terms  of  the  yield  stress  <ry,  the 
geometric  properties  of  the  section,  and  parameter  a. 

p  «•  2 atw<ry  (15a) 

m  —  m,  —  (15b) 


Equations  (15a)  and  (15b)  are  then  combined  to  eliminate  a,  resulting  in  the  following  para¬ 
bolic  equation  for  a  portion  of  the  yield  surface. 


■f 


UY 

4 iwm,  I  pr 


0 


(16) 


where 

A  section  width 

A  section  height 

A  section  area 

tf  flange  thtkness 
t9  web  thickness 


1 
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pp  plastic  force  with  sero  moment 
mf  plastic  moment  with  zero  force 


The  expressions  for  plastic  force  and  plastic  moment  are 

P,  -  Aw, 


jf-'/]]<v 


Case  2:  The  neutral  axis  in  flange.  (a  >  A/2  —  tj). 
By  a  similar  procedure  we  find  that 

P  *  Pr  ~  (*  ~ 


m 


-(*-4 


+  or 


b<rv 


which,  upon  combining  to  eliminate  a,  results  in 


m 

mm 


4  bmp 


,  _  jl\\  m  1+|  jl 

Pr  (  A  \Pr 


0 


(17) 

(18) 


(19a) 

(19b) 

(20) 


Together,  equations  (16)  and  (20)  define  the  yield  surface  for  positive  values  of  bending 
moment.  If  these  functions  are  plotted  using  the  dimensionless  variables 

-^-1  end  I— 


Pr 


m. 


the  function  is  symmetric  about  the  normalized  moment  axis.  When  the  sign  of  the  stresses  in 
the  assumed  distribution  is  reversed,  the  yield  functions  for  negative  moment  are  generated. 
Due  to  the  double  symmetry  of  the  cross  sections  under  consideration,  the  negative  moment 
yield  functions  are  simply  the  reflection  of  the  positive  moment  functions  about  the  normalized 
force  axis.  If  the  section  is  rectangular,  the  neutral  axis  is  always  in  the  web,  and  the  yield 
function  reduces  to  the  familiar  parabolic  function. 


When  the  positive  moment  yield  function  f*  is  plotted  together  with  the  negative 
moment  yield  function  f~,  there  is  a  discontinuity  in  the  slope  along  the  normalized  force  axis, 
forming  what  is  called  a  comer  condition.  To  obtain  a  continuous  variation  of  the  inelastic 
strain  rate  vector  as  the  corners  are  approached,  a  summation  of  the  two  yield  surfaces  will  be 
adopted  [6],  The  yield  functions  are  not  truncated  at  their  intersections,  but  extended 
indefinitely.  Three  types  of  region  are  thus  formed;  one  in  which  neither  function  is  active,  the 
second  in  which  only  one  function  is  active,  and  the  third  in  which  both  are  active.  The  resul¬ 
tant  strain  rate  vector  is  taken  as  the  vector  sum  of  the  contributions  from  all  active  yield  func¬ 
tions.  This  provides  a  continuous  variation  around  the  comer  region,  and  also  gives  a  vector 
parallel  to  the  force  axis  when  there  is  zero  bending  moment. 

It  is  necessary  to  discretize  the  constitutive  rate  equation  in  time,  and  to  linearize  it  with 
respect  to  the  primary  dependent  variables,  <  and  <r.  Following  the  procedure  employed  in  [2], 
the  rate  terms  are  replaced  with  an  Euler  backward  difference,  leading  to 

♦i+l  “  («J*i  ~  *.J  ~  C'jrrki  -  «r,|  +  yA<|<4(/>)>-^-  +  <*(/-)>  AC.|  (22) 


The  elastic  compliance  is 


The  vector  of  stresses  is 


C'  - 


—  0 
EA 

0  — 
El 


a  — 


The  vector  of  strains  for  small  displacements  is 


€  — 


Wax 

"\x* 


The  vector  of  strains  for  moderate  rotations  is 


u%x  + 


w  2 

W'X 


Taking  the  linear  part  of  eq.(22)  and  setting  it  equal  to  zero  defines  the  tangent  constitu¬ 
tive  law.  The  tangent  compliance,  C  r  .  will  contain  contributions  from  each  of  the  yield  func¬ 
tions  if  they  are  active.  For  the  specific  form,  see  [2]. 


4*&l]  ”  “  CrA<r «+>  “  0 


(23) 


Weighting  this  equation  with  a  virtual  stress  field  and  integrating  over  ft*  defines  the 
Galerkin  form. 


J^G'Jac  *+ 1  —  CrAai+i  +  “  0 


Gr  = 


bp  8m 


The  linearization  of  the  moderate  rotation  strain  vector, 

[Am,x  +  w,*Aw„j 

Ah,»«  |»4.j 


Ae;+|- 


(24a) 

(24b) 

(25) 


results  in  virtual  work  terms  which  are  conjugate  to  those  appearing  in  the  weak  form  of  equili¬ 
brium.  Hence,  the  symmetry  of  the  linearized  equations  is  preserved  for  moderate  rotations  as 
well  as  small  displacements. 

The  derivatives  of  /+  and  f~  required  in  eqs.(22)  and  (23)  are  conveniently  evaluated 
using  the  chain  rule. 


i  IL 

dp  pt  dip/pf) 

(26) 

ML.  i  */ 

dm  mr  dim/m,) 

(26b) 

To  keep  the  equations  dimensionally  consistent,  it  is  necessary  to  replace  fluidity  y  by  ypf 


4.  Finite  Elwwt  IrtwnMm 

In  order  that  the  integrals  in  the  Galerkin  forms  of  the  equilibrium  and  constitutive  equa¬ 
tions  may  be  evaluated,  it  is  necessary  to  specify  interpolation  functions  for  the  stress  resultants 
and  the  centroidal  axis  strains. 

The  displaced  shape  of  the  element  will  be  expressed  in  terms  of  C°  linear  shape  func¬ 
tions  for  axial  displacement  u(x),  and  C1  Hermite  polynomials  for  transverse  displacement 
v(x). 


N  ,(x)  0  0  Nj(x)  0  0  1  U} 

0  Hxix)  A/,(x)  0  Hi (x)  H4(x)  ]  UA 

Vi 


which  in  matrix  notation  becomes 


i(x)  -  N(x)U 


Strain  interpolation  relations  are  obtained  by  appropriately  differentiating  the  displacement 
shape  functions. 

For  small  displacements,  the  incremental  strain-displacement  operator  is  linear,  and 
results  in 


■(xlAU;*, 


-M 


For  moderate  rotations,  however,  the  operator  depends  upon  the  current  transverse  displace¬ 
ment  Aeld,  w(jt). 

(A  u,„  +  w.xAw,,! 

1  A.  =  |  (29) 

The  stresa  resultants  in  the  element  will  be  interpolated  using  piecewise  continuous  func¬ 
tions,  which  will  allow  for  discontinuities  in  the  stress  field  between  elements. 

Ip(x)  I  [l  0  0 1  p 

[ei(x)l~  a  .  JL  JL  mx 
L  L 


*  SOr)*:*,  (31) 

Substituting  these  interpolations  into  the  linearized  Oalerkin  forms  for  the  small  displace¬ 
ments  equations  leads  to 

8U'jf ■'0r)S0r)drAo(U1  -  8ujq„*,  -  J^B'(x)S(x)drw<U,|  (32) 


d*'[£s'(x)r(x)4rAI/J+i  -jfs'(x)CrS(x)drA«r 


S'Or)*^,*  (33) 
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E valuation  of  the  integrate  by  three  point  Gauss-Lobatto  numerical  quadrature  results  in 
the  following  mixed  system  of  equations 


0  T'  1  AU'J 
T  -qt]  a*;*. 


q,+1  - 

*!.+i 


(34) 


The  stress  interpolation  parameters  &  are  by  definition  piecewise  continuous  from  element  to 
element.  Hence,  the  mixed  system  may  be  statically  condensed  to  yield  a  generalized  displace¬ 
ment  model. 

T'QF*T  AU'+,  -  q„+i  -  T\r '+,  +  s  R'+1  (35a) 


Aa  -  -Q?  (♦'+,  -  TAU^+i)  (35b) 

In  the  moderate  rotation  case,  the  nonlinear  incremental  strain  displacement  matrix 
B(x,w)  must  be  used,  and  the  geometric  stiffness  K&  is  added  to  the  weak  form  of  equili¬ 
brium. 

(  Kc  T'(w) 

|T(w)  -Qr 

Statically  condensing, 

(kc  +  T'(h.)Q?‘T(h.))aU'+i  -  q,+i  -  T'(w)w'+1  +  T'(w)Qr'*J+i  5  R;+1  (37a) 

Adi +i  -  ~Qf‘  (♦;♦!  ~  T(w)AU'+1)  (37b) 


au;+1 

q,+i  -  T'UIo-J.m 

a*;+, 

(36) 


5.  Numerical  Results 

To  illustrate  the  selection  of  material  properties  for  the  viscoplastic  model,  two  sample 
problems  involving  the  axial  extension  of  a  one  element  bar  were  undertaken. 

In  the  first  example,  a  stress  relaxation  test,  the  bar  was  given  an  instantaneous  axial 
extension  and  the  axial  force  plotted  versus  time  for  various  values  of  the  fluidity  parameter  y. 
In  all  cases,  the  elastic  solution  was  taken  as  the  initial  condition.  The  results  of  this  test  are 
illustrated  in  Figure  3.  It  is  apparent  that  as  y  is  increased,  the  time  required  for  the  stresses  to 
relax  onto  the  yield  surface  decreases  in  an  exponential  manner.  Since  the  constitutive  equa¬ 
tion  is  of  necessity  discretized  in  time,  the  limiting  state  of  inviscid  plasticity  may  be  considered 
to  have  been  attained  when  the  time  required  for  relaxation  onto  the  yield  surface  is  less  than 
one  time  step  A  t.  The  value  of  y  corresponding  to  this  limiting  state  depends  upon  the 
material  properties  as  well  as  on  the  initial  stress  excursion  outside  of  the  yield  surface.  The 
user  is  warned  that  values  of  y  significantly  larger  than  this  limiting  value  may  result  in  physi¬ 
cally  improbable  solutions,  such  as  axial  contraction  of  an  element  when  subjected  to  a  tensile 
force  above  the  yield  load.  Hence,  it  is  imperative  that  proper  care  be  taken  in  specification  of 
the  material  properties. 

The  second  example  illustrates  the  effect  of  strain  rate  on  the  yield  stress.  The  end  of  the 
bar  was  extended  axially  at  constant  velocity,  and  axial  stress  plotted  versus  axial  strain  for  vari¬ 
ous  velocities.  The  value  of  y  was  held  constant  throughout.  Results  of  this  example  are  plot¬ 
ted  in  Figure  4.  When  the  strain  rate  is  increased,  the  effective  yield  stress  is  also  increased.  In 
the  limiting  case  of  vanishingly  small  strain  rate,  the  inviscid  plastic  response  is  recovered. 


The  inelastic  material  response  is  also  dependent  upon  the  value  of  the  exponent  used  in 
the  power  law  $(/).  A  similar  set  of  examples  should  also  be  undertaken  by  the  user  to  exam¬ 
ine  response  sensitivity  to  this  parameter.  It  has  been  found  by  the  authors  that  convergence 
difficulties  occur  if  the  exponent  iwwi  is  set  equal  to  1  whenever  bending  moment  gradients  are 
present  in  the  equilibrium  solution.  As  a  practical  consideration  therefore,  the  value  of  nnn 
should  always  be  set  equal  to  or  greater  than  2. 

To  illustrate  the  use  of  the  element  in  a  dynamic  analysis,  the  case  of  blast  loading  on  a 
simply  supported  beam  was  considered.  The  loading  was  idealized  as  a  uniformly  distributed 
triangular  pulse,  with  initial  amplitude  equal  to  2.5  times  the  static  collapse  load,  and  duration 
0.5  times  the  fundamental  period.  Due  to  the  symmetry  of  the  structure  and  loading,  only  one 
half  of  the  span  was  modeled,  using  ten  elements  of  equal  length.  Small  displacement  kinemat¬ 
ics  were  assumed.  The  value  of  y  was  set  to  model  inviscid  plasticity. 

Time  history  plots  were  made  of  the  midspan  deflection  and  bending  moment  (Figures  5 
A  6) .  Plots  were  also  made  of  the  distribution  of  bending  moment  along  the  half  span  at  every 
five  time  steps.  (Figures  7,8  A  9)  It  is  apparent  from  the  time  history  plots  that  the  response  is 
primarily  in  the  first  mode  of  the  structure.  The  period  of  constant  bending  moment  (yielding 
phase)  is  accompanied  by  a  parabolic  displacement  time  history.  After  sufficient  energy  is  dissi¬ 
pated  through  inelastic  deformation,  the  structure  responds  in  simple  harmonic  motion  at  the 
fundamental  frequency  about  a  residual  plastic  deformation.  The  ’noise'  which  appears  in  the 
moment  history  and  distributions  is  a  consequence  of  the  inertia  force  contribution  from 
accelerations  in  the  higher  modes. 

The  overall  response  of  the  beam  element  to  blast  loading  is  most  satisfactory,  and  com¬ 
pares  favorably  with  a  single  degree  of  freedom  approximation  given  in  (g)  . 
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Figure  1  •  Sign  Conventions 


Figure  2  •  Fully  Plastic  Strata  Distribution 
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Flgare  3  -  Axial  Force  Relaxation 
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Figare  4  -  Strain  Rate  Sensitivity  of  Yield  Stress 
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Flfut  5  -  Time  History  of  Midspan  Deflection 
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Figare  9  -  Bending  Moment  At  Time  0.03 5,  0.040  k  0.04S 
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Aypniii  I  -  C— psnen t  Matrices  Far  Ekanl  Efntltn 

The  beam  element  equations  are  given  in  general  form  by  eqs.(34)  and  (36),  and  contain 
three  typical  submatricea  which  result  from  the  weak  form  of  the  equilibrium  and  constitutive 
equations.  They  are  the  geometric  stiffness  Kc,  the  flexibility  Q,  and  the  transformation 
matrix  T.  Substitution  from  the  specified  interpolations  allows  these  matrices  to  be  integrated 
exactly,  resulting  in 

The  geometric  stiffness  (non  •  zero  terms  only): 

36  3£  —36  3Z, 

t  3  L  4  L2  -3  L  -L1 

*c  “  ~3oZ*  “36  -3 L  36  -31 
3 L  -L2  -3 L  4 L2 

The  flexibility  (elastic  case  only): 

0  0 

L  L 
3  El  6EI 
L  L 
6  El  3  El 

The  transformation  matrix  (small  displacements): 

-10  0  10  0 

T-  0  T  '10  i  0 

o  }  oo^-i 


The  transformation  matrix  (moderate  rotations): 

-la  b  1 
T(w)  -  0  y  -1  0 

0  1  0  0 


The  terms  a,  b,  c,  and  d  appearing  in  the  moderate  rotation  case  are  the  contraction  of  the 
geometric  stiffness  Kg  with  the  current  transverse  displacement  field  parameters  w. 
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Appeslix  II  -  Time  Integration  By  The  New  mark  Method 

The  incremental  dynamic  equations  of  equilibrium  for  a  discrete  system  are  written  in 
matrix  form  as 

MAA'*,  +  KrAU»+|  -  R.+i  -  MAa+| 

where  M  is  a  diagonal  lumped  mass  matrix(8]  ,  A  the  vector  of  nodal  accelerations,  KT  the 
tangent  stiffness  matrix,  U  the  nodal  displacement  vector  and  R  the  externally  applied  loads 
less  internal  resisting  forces,  all  evaluated  at  iteration  /  in  time  step  n+1. 

These  equations  are  coveniently  solved  using  the  Newmark  formulas  in  a  formulation 
having  displacement  as  primary  dependent  variable  [9].  Accordingly, 

U-i  -  U.  +  ArV,  +  -^(1  -  2fi) A, 

V.+1  -  V.  +  A/(l  -  y)A„ 

The  steps  for  a  typical  iteration  are  as  follows. 

Predictor  phase. 

u;+i  -6;+, 

v;+,  -  v„+1 

AJ+i  “  0 

Out  of  balance  forces. 

AB  “  RJ+i  —  MA^|  —  KfJJa+1 

Effective  stiffness  formulation. 

K  -  — L-M  +  Kr 

Ar^ 

K‘AU-  Aff 


Corrector  phase. 

u  ;:i  -  u  i+1  +  AU 
a;X1  -  (ll  ':l  -  Ua+iJ/fAr1^) 
V£l  -  V,+l  +  A/yA'tl 


Iterations  are  performed  until  a  satisfactory  convergence  tolerance  is  met.  The  time  step  is 
then  advanced  and  the  procedure  repeated. 
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Appendix  III  •  Beam  Element  Iayat  Data 

This  appendix  describes  the  input  data  necessary  to  use  the  plane  frame  viscoplastic  beam 
element  described  in  the  preceding  report.  General  input  data  required  for  use  of  program 
FEAP  is  outlined  in  chapter  24  of  (1],  and  will  not  be  covered  here.  However,  a  copy  of  the 
input  used  for  the  blast  load  analysis  example  appears  at  the  end  of  this  section. 

The  following  information  will  be  required  on  Card  2  of  the  general  input  data: 

Dimension  of  coordinate  space  •  2 
Degrees  of  freedom  per  node  -  3 
Nodes  per  element  -  2 
Added  degrees  of  freedom  •  0 

The  following  information  will  be  required  during  material  property  specification  following  the 
mate  macro  command  (3  cards): 

Card  1  -  format  (2il0) 

Field  1  -  material  set  number  as  specified  in  ekm  macro 
Field  2*1  (for  element  elmtOl) 

Card  2  -  format  (5x,a5,7fl0.0) 

Field  1  •  small  or  large  for  small  or  moderate  deformation  kinematics 
Field  2  •  Young's  modulus  E  (ksi) 

Field  3  *  section  width  b  (in) 

Field  4  -  section  depth  A  (in) 

Field  5  •  web  thickness  tm  (in) 

Field  6  -  flange  thickness  //  (in) 

Field  7  -  yield  stress  in  tension  cr y  (ksi) 

Field  8  •  fluidity  parameter  y  (1/sec) 

Card  3  -  (SflO.O) 

Field  1  •  exponent  nnn  with  minimum  value  of  2,  used  as  integer 
Field  2  -  mass  density  p  (kip-sec**2/in**4) 

Field  3  •  shift  of  yield  surface  origin  in  normalized  moment  direction 
Field  4  •  shift  of  yield  surface  origin  in  normalized  force  direction 
Field  5  -  tolerance  on  yield  criterion  /to I  (default  0.00001) 

The  shifts  of  the  yield  surface  origin  allow  the  user  to  move  the  yield  surface  along  the 
two  normalized  stress  axes.  The  tolerance  on  the  yield  criterion  allows  a  stress  state  within  ftol 
of  the  yield  surface  to  be  considered  elastic. 


Input  data  for  blast  loading  example  : 

feap  Sample  input  deck  for  simply  supported  beam 
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coor 

10 

1 

2 

3  2 

1 

1 

0.0 

0.0 

0.0 

11 

( blank  card) 

elem 

0 

50.0 

0.0 

0.0 

1 

( blank  card) 

boun 

1 

1 

2 

1 

1 

0 

1 

1 

0 

11 

( blank  card) 

fore 

0 

0 

0 

1 

2 

1 

0 

1.0 

0.0  0.0 

10 

0 

0 

1.0 

0.0  0.0 

11 

( blank  card) 

mate 

1 

0 

1 

0 

0.5 

0.0  0.0 

smalt  30000.0 

2.0 

10.0 

5.0  40. 

2.0  0.00001 

( blank  card) 

0.0 

0.0 

0.0 

